function [notchfilt,chebyfilt,sgofilt] = ecgfilt(data)
    %% Filtering
    %% Load data
    %Initializing variables
    Samp_freq = 553;

    %% Setting cut-off frequencies
    %Hi_cutoff = 2*0.5/Samp_freq; %High-pass filtfilt
    %Lo_cutoff = 2*99/Samp_freq; %Low-pass filtfilt
    Stop_1 = 2*49.005/Samp_freq; %Bandstop cutoff freqencies
    Stop_2 = 2*50.005/Samp_freq;
    %changing from 16bits to double for filtfilting analysis
    double_data = double(data);

    %% Bandstop filtfilt
    [F,E] = butter(4,[Stop_1,Stop_2],'stop');
    x_stop = filtfilt(F,E,double_data);
    notchfilt = x_stop;
    
    %% Chebyshev ripple lowpass filtfilt
    %Rp: pass band ripple in dB
    %Rs: stop band ripple attenuation in dB
    %Wp: pass band
    %Ws: stop band in ripples
    Wp1 = 96 / (Samp_freq/2);
    Ws1 = 100 / (Samp_freq/2);
    Rp1 = 5;
    Rs1 = 30;
    %Domain Transformation
    [n,Wp1] = cheb1ord(Wp1,Ws1,Rp1,Rs1);
    %Domain re-transform
    [D1,C1] = cheby2(n,Rp1,Wp1);
    x_cheby1 = filtfilt(D1,C1,x_stop);
    chebyfilt = x_cheby1;
    %% Savitzky-Golay filtfilting
    %order: Wave preciseness
    %windows: Frame preciseness
    order = 3;
    windows = 31;
    x_sgo = sgolayfilt(x_cheby1, order, windows);
    sgofilt = x_sgo;
